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Abstract: The Heston model stands out from the class of stochastic volatility (SV) models mainly for 
two reasons. Firstly, the process for the volatility is non-negative and mean-reverting, which is what we 
observe in the markets. Secondly, there exists a fast and easily implemented semi-analytical solution for 
European options. In this article we adapt the original work of Heston (1993) to a foreign exchange (FX) 
setting. We discuss the computational aspects of using the semi-analytical formulas, performing Monte 
Carlo simulations, checking the Feller condition, and option pricing with FFT. In an empirical study we 
show that the smile of vanilla options can be reproduced by suitably calibrating three out of five model 
parameters. 
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1.1 Introduction 

The universal benchmark for option pricing is flawed. The Black-Scholes formula is based on the 
assumption of a geometric Brownian motion (GBM) dynamics with constant volatility. Yet, the 
model-implied volatilities for different strikes and maturities of options are not constant and tend 
to be smile shaped (or in some markets skewed). Over the last three decades researchers have tried 
to find extensions of the model in order to explain this empirical fact. 

As suggested already by Merton (1973), the volatility can be made a deterministic function of time. 
While this approach explains the different implied volatility levels for different times of maturity, it 
still does not explain the smile shape for different strikes. Dupire (1994), Derman and Kani (1994), 
and Rubinstein (1994) came up with the idea of allowing not only time, but also state dependence 
of the volatility coefficient, for a concise review see e.g. Fengler (2005). This local (deterministic) 
volatility approach yields a complete market model. It lets the local volatility surface to be fitted, 
but it cannot explain the persistent smile shape which does not vanish as time passes. Moreover, 
exotics cannot be satisfactorily priced in this model. 

The next step beyond the local volatility approach was to allow the volatility to be driven by a 
stochastic process. The pioneering work of Hull and White (1987), Stein and Stein (1991), and 
Heston (1993) led to the development of stochastic volatility (SV) models, for reviews see Fouque, 
Papanicolaou, and Sircar (2000) and Gatheral (2006). These are multi- factor models with one 
of the factors being responsible for the dynamics of the volatility coefficient. Different driving 
mechanisms for the volatility process have been proposed, including GBM and mean-reverting 
Ornstein-Uhlenbeck type processes. 
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The Heston model stands out from this class mainly for two reasons. Firstly, the process for the 
volatility is non-negative and mean-reverting, which is what we observe in the markets. Secondly, 
there exists a fast and easily implemented semi-analytical solution for European options. This 
computational efficiency becomes critical when calibrating the model to market prices and is the 
greatest advantage of the model over other (potentially more realistic) SV models. Its popularity 
also stems from the fact that it was one of the first models able to explain the smile and simulta- 
neously allow a front-office implementation and a valuation of many exotics with values closer to 
the market than the Black-Scholcs model. Finally, given that all SV models generate roughly the 
same shape of implied volatility surface and have roughly the same implications for the valuation 
of exotic derivatives (Gatheral, 2006), focusing on the Heston model is not a limitation, rather a 
good staring point. 

This chapter is structured as follows. In Section [TT2"1 we define the model and discuss its properties, 
including marginal distributions and tail behavior. Next, in Section [T751 wc adapt the original work 
of Heston (1993) to a foreign exchange (FX) setting. We do this because the model is particularly 
useful in explaining the volatility smile found in FX markets; in equity markets the typical volatility 
structure is a strongly asymmetric skew (also called a smirk or grimace). In Section [L4l we show 
that the smile of vanilla options can be reproduced by suitably calibrating the model parameters. 
Finally, in Section [T31 we briefly discuss the alternatives to the Heston model. 



1.2 The model 

Following Heston (1993) consider a stochastic volatility model with GBM-like dynamics for the 
spot price: 




and a non-constant instantaneous variance Vt driven by a mean-reverting square root (or CIR) 
process: 

dv t = n{9 - v t ) dt + ojv~ t dwf> '. (1.2) 

The stochastic increments of the two processes are correlated with parameter p, i.e. dW^ dW^ = 
pdt. The remaining parameters - fi, 9, k, and a - can be interpreted as the drift, the long-run 
variance, the rate of mean reversion to the long-run variance, and the volatility of variance (often 
called the vol of vol), respectively. Sample paths and volatilities of the GBM and the Heston spot 
price process arc plotted in Figure fTTTI 

By setting x t = log (5* /So) — M*: we can express the Heston model (|1.I|) - (|I.2|) in terms of the 
centered (log-)return xt and Vt- The process is then characterized by the transition probability 
Pt(x,v \vo) to have (log-)return x and variance v at time t given the initial return x — and 
variance vq at time t = 0. The time evolution of Pt(x,v | vo) is governed by the following Fokker- 
Planck (or forward Kolmogorov) equation: 



dt 



P = 



^{{v-e)p} 



id_ 

2dx 
I d 2 



(vP) 



d < m 1 9 , ox d < on 

+ P ^V (VP)+ 2^ VP) + T^ VP) - 



(1.3) 



Solving this equation yields the following semi-analytical formula for the density of centered returns 
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Figure 1.1: Sample trajectories (left panel) and volatilities (right panel) of the GBM and the Heston 
spot price process (|1.1[) obtained with the same set of random numbers. 
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x, given a time lag t of the price changes (Dragulescu and Yakovenko, 2002): 

P t (x) = — / e^+^d^ (1.4) 

with 

F t (0 = g * - log (cosh f + sinh f ) , 

7 = « + and 51 = \f^f 2 + er 2 (£ 2 — i£). 

Somewhat surprisingly, the introduction of SV does not change the properties of the spot price 
process in a way that could be noticed just by a visual inspection of its realizations, see Figure 
ll.ll where sample paths of a GBM and the spot process (jl.ip in the Heston model are plotted. To 
make the comparison more objective both trajectories were obtained with the same set of random 
numbers. 

In both cases the initial spot rate So = 4 and the domestic and foreign interest rates are = 5% 
and rf = 3%, respectively, yielding a drift of fj, = r<i — rf = 2%. The volatility in the GBM is 
constant ^Jv[ = %/ 4% = 20%, while in the Heston model it is driven by the mean- reverting process 
(|1.2[) with the initial variance vq = 4%, the long-run variance 9 = 4%, the speed of mean reversion 
ft = 2, and the vol of vol a = 30%. The correlation is set to p = —0.05. 

A closer inspection of the Heston model does, however, reveal some important differences with 
respect to GBM. For instance, the probability density functions (pdfs) of (log-)returns have heavier 
tails - exponential compared to Gaussian, see Figure 11.21 In this respect they are similar to 
hyperbolic distributions (Weron, 2004, see also Chapter 1), i.e. in the log-linear scale they resemble 
hyperbolas, rather than parabolas of the Gaussian pdfs. 
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Figure 1.2: The marginal pdfs in the Black-Scholes (GBM) and Heston models for the same set 
of parameters as in Figure H~T1 (left panel). The tails of the Heston marginal pdfs are 
exponential, which is clearly visible in the log-linear scale (right panel). 
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1.3 Option pricing 

Consider the value function of a general contingent claim U(t,v,S) paying g(S) = U(T,v,S) at 
time T. We want to replicate it with a self-financing portfolio. Due to the fact that in the Heston 
model we have two sources of uncertainty (the Wiener processes and W^) the portfolio 

must include the possibility to trade in the money market, the underlying and another derivative 
security with value function V(t, v, S). 

We start with an initial wealth Xq which evolves according to: 

dX = A dS + T dV + r d (X - TV) dt - (r d - r f )AS dt, (1.5) 

where A is the number of units of the underlying held at time t and T is the number of derivative 
securities V held at time t. 

The goal is to find A and T such that X t = U(t,v t ,St) for all t £ [0, T]. The standard approach 
to achieve this is to compare the differentials of U and X obtained via Ito's formula. After some 
algebra we arrive at the partial differential equation (PDE) which U must satisfy (for details on 
the derivation in the foreign exchange setting see Hakala and Wystup, 2002): 

1 a2 d 2 U a d 2 U 1 2 d 2 U , .BU 

2 vS as* + pavS dsd-v + r v ^ + {rd - rf)s ds + 

+{ K (9-v)-X(t,v,S)}^-r d U + ^ = 0. (1.6) 



The term X(t, v, S) is called the market price of volatility risk. Heston (1993) assumed it to be linear 
in the instantaneous variance Vt, i.e. \(t, v, S) = Xvt, in order to retain the form of the equation 
under the transformation from the statistical (or risky) measure to the risk-neutral measure. 
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1.3.1 European vanilla FX option prices and Greeks 

We can solve (|1.6[) by specifying appropriate boundary conditions. For a European vanilla FX 
option these are: 

U(T,v,S) = max{(f>(S - K),0}, (1.7) 
U(t,v,0) = *-l*Ke- r < T , (1.8) 

oo) = l±*c^, (1.9) 

r d U(t,0,S) = (r d -r / )^||(t,0,5) + 

+ K0^(i,O,S) + ^(t,O,S), (1.10) 
rrr* ctt JSe"^, for = +1, 

f7 (t,oo,S) = < 1.11 
I Ke rdT , for = -1, 

where <f> = ±1 for call and put options, respectively. The strike is in units of the domestic 
currency and r = T — t is the time to maturity (i.e. T is the expiration time in years and t is the 
current time). 

Heston (1993) solved the PDE analytically using the method of characteristic functions. For 
European vanilla FX options the price is given by: 

h(r) = HestonVanilla(K, 9, er, p, A, rf, vt, St, K, r, cb) 

= ^{e-^StP+i^-Ke-^P.^)}, (1.12) 

where Ui t 2 = ±5, 61 = K + A — op, b 2 = « + A, and 



^/ (poipi — 6j) 2 — o 2 (2ujtpi — ip 2 ), (1-13) 

ft = ^-^ + ^ , (1.14) 

bj — paipt — clj 

Cj(T,ip) = {r d -r f )ipiT+ (1-15) 
+ — \ (bj - paipi + d 3 )T - 2 lor ' 



6j — per (pi + dj ( 1 — e djT 
l-^-e^ 7 



^) = J T T M t^). (1-16) 



fj(x,v t ,T,<p) = exp{Cj(T,<p) + Dj(r,(p)v t + itpx}, (1.17) 
1 , 1 /"~»/e~ <w /j(*. «*>*■,¥>) 



Pjfowt.r.y) = dy. (1-18) 

Note that the functions Pj are the cumulative distribution functions (in the variable y = log K) of 
the log-spot price after time r = T — t starting at x = log St for some drift p. Finally: 

P+{4>) = ^^+<t>Pi(x,v U T,y), (1.19) 

P-{4>) = ^-^ + cbP 2 (x,v t ,T,y). (1.20) 
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The Greeks can be evaluated by taking the appropriate derivatives or by exploiting homogeneity 
properties of financial markets (Reiss and Wystup, 2001). In the Heston model the spot delta and 
the so-called dual delta are given by: 



dh{t) 



r P+{<t>) and 



dh(t) 



r P-(4>), 



dSt OK 
respectively Gamma, which measures the sensitivity of delta to the underlying has the form 

„ dA e~ r f T 



(1.21) 



dS t 



Si 



-pi(logSt,v t ,T,logK), 



where 



Pj(x,v,T,y) = - 3?{e %vv fj{x,v,T,<p)} dip, j = 1,2, 



(1.22) 
(1.23) 



are the densities corresponding to the cumulative distribution functions Pj (|1.18l) . 

The time sensitivity parameter theta = dh(t)/dt can be easily computed from (jl.6p . while the 
formulas for rho are the following: 



dh(t) 
dr d 



(f>Ke~ riT rP-{4>), 



dh(t) 
dr 



f 



-cj>Ste- r f r TP + {(j>). 



(1.24) 



Note, that in the foreign exchange setting there are two r/io's - one is a derivative of the option 
price with respect to the domestic interest rate and the other is a derivative with respect to the 
foreign interest rate. 

The notions of vega and volga usually refer to the first and second derivative with respect to 
volatility. In the Heston model we use them for the first and second derivative with respect to the 
initial variance: 



dh(t) 
dv t 



8 2 h(t) 

dv? 



e-^St—PtilogSuVt^logK) - 
ov t 

- Ke rdT A P 2 (log S t ,v u r, log K) , 

OVt 

e-^St^Px (log S u v t , r, log K) - 
ovt 



where 



Pj(x,v t ,T, y) 



^Pj(x,v U T, y) 



Ke~ 



'-^P 2 (logS t ,v t ,T, log K) , 



D(T,<p)e-*nfj(x,v t ,T,(p) 



D 2 (T,ip)e-^yf,(x,v t ,T,ip) 



ip 



d<p, 
d(p. 



(1.25) 

(1.26) 

(1.27) 
(1.28) 



1.3.2 Computational issues 

Heston's solution is semi-analytical. Formulas (|1.19ll.2U|l require to integrate functions fj, which 
are typically of oscillatory nature. Hakala and Wystup (2002) propose to perform the integration 
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in (|1.18p with the Gauss-Laguerre quadrature using 100 for oo and 100 abscissas. Jackel and Kahl 
(2005) suggest using the Gauss-Lobatto quadrature (e.g. Matlab's quadl.m function) and transform 
the original integral boundaries [0, +oo) to the finite interval [0, 1]. 

As a number of authors have recently reported (Albrecher et al., 2006; Gathcral, 2006; Jackel 
and Kahl, 2005), the real problem starts when the functions fj arc evaluated as part of the 
quadrature scheme. In particular, the calculation of the complex logarithm in cqn. (|1.15p is prone 
to numerical instabilities. It turns out that taking the principal value of the logarithm causes Cj 
to jump discontinuously each time the imaginary part of the argument of the logarithm crosses 
the negative real axis. One solution is to keep track of the winding number in the integration 
(|1.18|) . but is difficult to implement because standard numerical integration routines cannot be 
used. Jackel and Kahl (2005) provide a practical solution to this problem. 

A different approach was taken by Albrecher et al. (2006), see also Gatheral (2006). They proposed 
to make a simple transformation when computing the characteristic function and proved that 
the numerical stability of the resulting formulas is guaranteed under the full dimensional and 
unrestricted parameter space. Namely, their idea is to switch from gj in (|1.14j) to 

gj bj — potpi + dj 
which leads to new formulas for Cj and Dj-. 

Cj(T,<p) = {r d -r f )<piT + (1.30) 

+ ^2 | W - P ai P l - d 3) T - 2 lo S ^ — 

Dj(r,<p) = 3 ^ I °- ), 1.31 

c \ 1 — gj e a i T j 

in p,17p . Note, that the only differences between eqns. (|1.14j) - (|1.16[) and (|1.29l) - (|1.31l) . respectively, 
are the flipped minus and plus signs in front of the dj's. 

The mispricings resulting from using (|1.14[) - (jl.l6p are not that obvious to notice. In fact, if we 
price and backtest on short or middle term maturities only, we might not detect the problem at 
all. However, the deviations can become extreme for longer maturities (typically above 3-5 years; 
the exact threshold is parameter dependent, see Albrecher et al., 2006). 

Apart from the above semi-analytical solution for vanilla options, alternative approaches can be 
utilized. These include the general Fast Fourier Transform (FFT) approach of Carr and Madan 
(1999), finite differences, finite element methods and Monte Carlo simulations. The FFT-based 
pricing technique is discussed in Section 11.3.41 As for the other methods, finite differences must 
be used with care since high precision is required to invert scarce matrices. The Crank-Nicholson, 
ADI (Alternate Direction Implicit), and Hopscotch schemes can be used, however, ADI is not 
suitable to handle nonzero correlation. Boundary conditions must be set appropriately, for details 
see Kluge (2002). On the other hand, finite element methods can be applied to price both the 
vanillas and exotics, as explained for example in Apel, Winkler, and Wystup (2002). 

Finally, the Monte Carlo approach requires attention as the simple Euler discretization of the 
CIR process (jl.2p may give rise to a negative variance. To deal with this problem, practi- 
tioners generally either adopt (i) the absorbing v t — max(w t ,0), or (ii) the reflecting principle 
v t = max(i!i, — vt). More elegant, but computationally more demanding solutions include sam- 
pling from the exact transition law (Glasserman, 2004) or application of the quadratic-exponential 
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(QE) scheme (Andersen, 2008). For a recent survey see the latter reference, where several new 
algorithms for time-discretization and Monte Carlo simulation of Hcston-type SV models are in- 
troduced and compared in a thorough simulation study. Both the absorbing/reflecting patches and 
the QE scheme are implemented in the simHeston.m function used to plot Figure ITTT1 



1.3.3 Behavior of the variance process and the Feller condition 

The CIR process for the variance, as defined by (|1.2p . always remains non- negative. This is 
an important property which, for instance, is not satisfied by the Ornstein-Uhlenbeck process. 
However, ideally one would like a variance process which is strictly positive, because otherwise the 
spot price process degenerates to a deterministic function for the time the variance stays at zero. 
As it turns out, the CIR process remains strictly positive under the condition that 

a := 4 4 > 2, (1.32) 

which is often referred to as the Feller condition. We call a the dimensionality of the corresponding 
Bessel process (see below). If the condition is not satisfied, i.e. for < a < 2, the CIR process will 
visit recurrently but will not stay at zero, i.e. the 0-boundary is strongly reflecting. 

Unfortunately, when calibrating the Heston model to market option prices it is not uncommon 
for the parameters to violate the Feller condition (|1.32p . This is not a complete disaster, as the 
variance process can only hit zero for an infinitesimally small amount of time, but it is still worrying 
as very low levels of volatility (e.g. say below 1%) are repeatedly reached for short amounts of time 
and that is not something observed in the market. 

Besides being important from a modeling point of view, the Feller condition also plays a role 
in computational accuracy. For Monte Carlo simulations special care has to be taken so that the 
simulated paths do not go below zero if (|1.32p is not satisfied. On the PDE side, the Feller condition 
determines whether the zero-variance boundary is in- or out-flowing, that is to say whether the 
convection vector at the boundary points in- or outwards. To see this, we need to write the log-spot 
transformed Heston PDE in convection-diffusion form 

d 

— U = div(Agrad U) - dw{Ub) + f, (1.33) 



and obtain 



which is out-flowing at the v — boundary if 

icr 2 - K 6 < 0, (1.35) 
in which case the Feller condition (|1.32j) is satisfied. 

Having introduced and discussed the importance of this condition, we now give an overview of 
how it is derived; we refer the interested reader to Chapter 6.3 in Jeanblanc, Yor, and Chesney 
(2009) for a more thorough treatment. The main idea is to relate the CIR process to the family of 
squared Bessel processes which have well known properties. We call X t an a-dimensional squared 
Bessel process and denote it by BES 2 (a) if it follows the dynamics of 



dX t =adt + 2JX+ dW t , (1.36) 
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with := max{0,X t }. This definition makes sense for any real valued a. However, in the case 
of integer valued a we have an interesting interpretation: a BES 2 (a) process X t , with Xq = 0, 
follows the same dynamics as the squared distance to the origin of an a-dimensional Brownian 
motion, i.e. 

a 

X t = Y / B i t l \ (1.37) 

i=l 

with being independent Brownian motions. From this we can conclude that for a = 2, X t will 
never reach zero and, using the stochastic comparison theorem (Rogers and Williams, 2000), this 
property remains true for any a > 2. Similarly for < a < 1 the value zero will be repeatedly hit 
(for a = 1 this happens as often as a one-dimensional Brownian motion crosses zero). 

The stochastic comparison theorem also immediately tells us that BES 2 (a) processes are non- 
negative for non-negative a: for a — we get the trivial solution X t = and increasing the drift 
term to a > cannot make the paths any smaller. For a proof of the above statements we refer 
to Chapter V.48 in Rogers and Williams (2000), from which we also state the following additional 
properties. If X t is an a-dimensional squared Bessel process BES 2 (a) and Xq > then it is always 
non-negative and: 

1. for < a < 2 the process hits zero and this is recurrent, but the time spent at zero is zero, 

2. for a = 2 the process is strictly positive but gets arbitrarily close to and oo, 

3. for a > 2 the process is strictly positive and tends to infinity as time approaches infinity. 

To translate these properties to the class of CIR processes, we only need to apply the following 
space-time transformation to the squared Bessel process. Define dY t := c~ Kt Xp^t Then Y t 
follows the dynamics of 

dY t = K(a/3 - Y t ) dt + 2^/n/3Y t dW t , (1.38) 
which is the same as the dynamics of the CIR process (|1.2|) if we set j3 = and a = 4 = ■ 



1.3.4 Option pricing with FFT 

In this section, we briefly describe the numerical option pricing approach of Carr and Madan (1999), 
which utilizes the characteristic function (cf) of the underlying instrument's price process. The 
basic idea of the method is to compute the price by Fourier inversion given an analytic expression 
for the cf of the option price. 

The rationale for using this approach is twofold. Firstly, the algorithm offers a speed advantage, 
including the possibility to calculate prices for a whole range of strikes in a single run. Secondly, 
the cf of the log-price is known and has a simple form for many models, while the pdf is often 
either unknown in closed-form or complicated from the numerical point of view. For instance, for 
the Heston model the cf takes the form (Heston, 1993; Jackel and Kahl, 2005): 

E{exp(^logS* T )} = f 2 (x, v t , r, <p), (1.39) 

where f 2 is given by (|1.17[) . 

Let us now derive the formula for the price of a European vanilla call option. Derivation of the put 
option price follows the same lines, for details see Lee (2004) and Schmelzlc (2010). Alternatively 
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we can use the call-put parity for European vanilla FX options (see e.g. Wystup, 2006). Let /i c (t; k) 
denote the price of the call option maturing in r = T — t years with a strike of K — exp(fc): 



h c ( T ;k)=j e- rl (e s - e k )q T (s)ds, (1.40) 

where qt is the risk-neutral density of st = log St- The function h c (r; k) is not square-intcgrablc 
(see e.g. Rudin, 1991) because it converges to Sq for k —> — oo. Hence, consider a modified function 
H c (r;k) = e ak h c (r; k), which is squarc-intcgrable for a suitable constant a > 0. A sufficient 
condition is given by: 

E{(5 T ) Q+1 }<oo, (1.41) 

which is equivalent to V't(O), i.e. the Fourier transform of H (r; fc), see (|1.42[) below, being finite. 
In an empirical study Schoutens, Simons, and Tistaert (2004) found that a = 0.75 leads to stable 
algorithms, i.e. the prices are well replicated for many model parameters. This value also fulfills 
condition (|1.41j) for the Heston model (Borak, Detlefsen, and Hardle, 2005). Note, that for put 
options the condition is different: it is sufficient to choose a > such that E{(Sx)~ a } < oo 
(Lee, 2004). 

Now, compute the Fourier transform of H c (t; k): 

/>OG 

ifr(v) = I e lvk H c { T -k)dk 

p OO 

„ufe / ak —rT / s k\ ( „\ j „jj 

e / e e (e — e )qT\s)dsdk 



e~ rl Qt(s) / | e a * +s ~ e \ a + l >* | e vv «dkds 

f — oo 

Ja+l+iv)s (a.+l+iv)s 
-rT , / n I c 



oo 



e—grOO — 77— 

|_ a + iu a + l + zf 

^/ 2 { v -(a + l)i} 



(1.42) 



a 2 + a — v 2 + i(2a + l)u ' 
where / 2 is the cf of qx, see (JT3HJ). We get the option price in terms of ipT using Fourier inversion: 

h c {T . k) = cx P(~ afc ) r e - ivk i>(v)dv. (1.43) 
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This integral can be numerically approximated as (Carr and Madan, 1999): 

h c (r;k u ) «^Ve-*(' , - 1 )M)e««^( ( .)!{3|(-l/-{ J , 1 }, (1.44) 

7T ' — ' O 

where fc u = i{— 7r + ^-(u — 1)}, u = 1, . . . , JV", 17 > is the distance between the points of the 
integration grid, Vj = r)(j — 1), j = 1, . . . , N, and <5 is the Dirac function. 



As can be seen in Figure fL3l the differences between the FFT method and the (semi-)analytical 
formula (|1.12j) arc relatively small. The differences result form the fact that the former method 
yields 'exact' values only on the grid k u . In order to preserve the speed of the FFT-based algorithm 
we use linear interpolation between the grid points. This approach, however, slightly overestimates 
the true prices since the option price is a convex function of the strike (Borak, Detlefsen, and 
Hardle, 2005). It can be clearly seen that near the grid points the prices obtained by both methods 
coincide, while between the grid points the FFT-bascd algorithm generates higher prices than the 
analytical solution, see the right panels in Figure 11.31 
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Figure 1.3: European call (top left) and put (bottom left) FX option prices obtained using the FFT 
method and the 'Analytical' formula (|1.12p for a sample set of parameters. Right panels: 
The percentage differences between the prices: (FFT — Analytical)/ Analytical x 100%. 
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1.4 Calibration 

1.4.1 Qualitative effects of changing the parameters 

Before calibrating the model to market data we will show how changing the input parameters 
affects the shape of the fitted smile curve. This analysis will help in reducing the dimensionality 
of the problem. In all plots of this subsection the solid blue curve with x's is the smile obtained 
for v = 0.01, a = 0.2, re = 1.5, 9 = 0.015, and p = 0.05. 

First, take a look at the initial variance (top left panel in Figure [L^|l . Apparently, changing vq 
allows for adjustments in the height of the smile curve. On the other hand, the volatility of variance 
(vol of vol) has a different impact on the smile. Increasing a increases the convexity of the fit, see 
the top right panel in Figure 11.41 In the limiting case, setting a equal to zero would produce a 
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Figure 1.4: The effects of changing the model parameters on the shape of the smile: initial variance 
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deterministic process for the variance and, hence, a volatility which does not admit any smile (a 
constant curve). 

The effects of changing the long-run variance 9 are similar to those observed by changing the initial 
variance, compare the left panels in Figure 11.41 It seems promising to choose the initial variance 
a priori, e.g. set ^/vq = implied at-the-money (ATM) volatility, and only let the long-run variance 
9 vary. In particular, a different initial variance for different maturities would be inconsistent. 

Changing the mean reversion k affects the ATM part more than the extreme wings of the smile 
curve. The low/high deltas (A) remain almost unchanged while increasing the mean reversion 
lifts the center, see the bottom right panel in Figure 11.41 Moreover, the influence of k is often 
compensated by a stronger vol of vol a. This suggests fixing mean reversion (at some level, say 
K = 1.5) and only calibrating the remaining three parameters. If the obtained parameters do not 
satisfy the Feller condition (|1.32[) . it might be worthwhile to fix k at a different level, say ft = 3, 
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recalibrate the remaining parameters and check if the new estimates fulfill the condition and lead 
to a more realistic variance process. 

Finally, let us look at the influence of correlation. The uncorrelated case produces a fit that looks 
like a symmetric smile curve centered at-the-money, see Figure 11.51 However, it is not exactly 
symmetric. Changing p changes the degree of symmetry. In particular, positive correlation makes 
calls more expensive, negative correlation makes puts more expensive. Note, that the model 
yields a volatility skew, a typically observed volatility structure in equity markets, only when the 
correlation is set to a very high or low value 

1.4.2 The calibration scheme 

Calibration of SV models can be done in two conceptually different ways. One is to look at the 
time series of historical data. Estimation methods such as Generalized, Simulated, and Efficient 
Methods of Moments (respectively GMM, SMM, and EMM), as well as Bayesian MCMC have 
been extensively applied. See Broto and Ruiz (2004) for a review. In the Heston model we could 
also try to fit empirical distributions of returns to the marginal distributions specified in (|1.4[) via 
a minimization scheme. Unfortunately, all historical approaches have one common flaw: they do 
not allow for the estimation of the market price of volatility risk X(t,v,S). Observing only the 
underlying spot price and estimating SV models with this information will not yield correct prices 
for the derivatives. 

This leads us to the second estimation approach: instead of using the spot data we calibrate the 
model to the volatility smile (i.e. prices of vanilla options). In this case we do not need to worry 
about estimating the market price of volatility risk as it is already embedded in the market smile. 
This means that we can set A = by default and just determine the remaining parameters. 

As a preliminary step, we have to retrieve the strikes since the smile in FX markets is specified as 
a function of delta. Comparing the Black-Scholcs type formulas (in the FX market setting we have 
to use the Garman and Kohlhagen (1983) specification) for delta and the option premium yields 
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the relation for the strikes Ki. From a computational point of view this stage requires only an 
inversion of the Gaussian distribution function. Next, based on the findings of Section [1.4. 11 we fix 
two parameters (initial variance vq and mean reversion k) and fit the remaining three: volatility of 
variance a, long-run variance 6, and correlation p for a fixed time to maturity and a given vector 
of market Black-Scholes implied volatilities {<7i}f =1 for a given set of delta pillars {Ai}™ =1 . 

After fitting the parameters we compute the option prices in the Heston model using (|1.12p and 
retrieve the corresponding Black-Scholes model implied volatilities {o"i}™ =1 via a standard bisection 
method (function fzero.m in Matlab uses a combination of bisection, secant, and inverse quadratic 
interpolation methods). The next step is to define an objective function, which we choose to be 
the Sum of Squared Errors (SSE): 



We compare the volatilities because they are all of similar magnitude, in contrast to the prices 
which can range a few orders of magnitude for it-the-money (ITM) vs. out-of-the-money (OTM) 
options. In addition, one could introduce weights for all the summands to favor ATM or OTM fits. 
Finally we minimize over this objective function using a simplex search routine (fminsearch.m in 
Matlab) to find the optimal set of parameters. 

1.4.3 Sample calibration results 

We are now ready to calibrate the Heston model to market data. First, we take the EUR/USD 
volatility surface on July 1, 2004 and fit the parameters in the Heston model according to the 
calibration scheme discussed earlier. The results are shown in Figures 11.6111.71 Note, that the fit 
is very good for intermediate and long maturities (three months and more). Unfortunately, the 
Heston model does not perform satisfactorily for short maturities (see Section H. 5. 2 1 for a discussion 
of alternative approaches). Comparing with the fits in Weron and Wystup (2005) for the same 
data, the long maturity (2Y) fit is better. This is due to a more efficient optimization routine 
(Matlab 7.2 vs. XploRe 4.7) and utilization of the transformed formulas (|1.29[) - (|1.31|) instead of 
the original ones (|1.14l) - (|1.16l) . 

Now we take a look at more recent data and calibrate the Heston model to the EUR/USD volatility 
surface on July 22, 2010. The results are shown in Figures fl. 81 11.91 Again the fit is very good 
for intermediate and long maturities, but unsatisfactory for maturities under three months. The 
term structure of the vol of vol and correlation visualizes the problem with fitting the smile in 
the short term for both datasets. The calibrated vol of vol is very low for the 1W smiles, then 
jumps to a higher level. The correlation behaves similarly, for the 1W smiles it is much lower than 
for the remaining maturities. In 2004 it additionally changes sign as the skew changes between 
short and longer-term tenors. Also note, that the more skewed smiles in 2010 require much higher 
(anti-)correlation (—0.4 < p < —0.3 for r > 1M) to obtain a decent fit than the more symmetric 
smiles in 2004 (-0.01 < p < 0.05 for t > 1M). 

As these examples show, the Heston model can be successfully applied to modeling the volatility 
smile of vanilla FX options in the mid- to longer-term. There are essentially three parameters to 
fit, namely the long-run variance (#), which corresponds to the ATM level of the market smile, 
the vol of vol (tr), which corresponds to the convexity of the smile (in the market often quoted as 
butterflies), and the correlation (p), which corresponds to the skew of the smile (quoted as risk 




(1.45) 
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Figure 1.6: The EUR/USD market smile on July 1, 2004 and the fit obtained with the Heston 
model for r = 1 week (top left), 1 month (top right), 3 months (bottom left), and 6 
months (bottom right). 
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reversals). It is this direct link of the model parameters to the market that makes the Heston 
model so attractive to practitioners. 

The key application of the model is to calibrate it to vanilla options and afterward use it for pricing 
exotics (like one-touch options, see Wystup, 2003) in either a finite difference grid or a Monte Carlo 
simulation. Surprisingly, the results often coincide with the traders' rule of thumb pricing method 
(Wystup, 2006). This might also simply mean that a lot of traders use the same model. After 
all, it is a matter of belief which model reflects reality the best. Recent ideas are to take prices of 
one-touch options along with prices of vanilla options from the market and use this common input 
to calibrate the Heston model. 
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Figure 1.7: The EUR/USD market smile on July 1, 2004 and the fit obtained with the Heston 
model for r = 1 year (top left) and 2 years (top right). The term structure of the vol 
of vol and correlation visualizes the problem with fitting the smile in the short term 
(bottom panels). 
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1.5 Beyond the Heston model 
1.5.1 Time-dependent parameters 

As we have seen in Section 11.4.31 performing calibrations for different time slices of the volatility 
matrix produces different values of the parameters. This suggests a term structure of some param- 
eters in the Heston model. Therefore, we need to generalize the CIR process (|1.2p to the case of 
time-dependent parameters, i.e. we consider the process: 



dv t = n(t){0(t) - v t } dt + a(t)y/v t dW t , 



(1.46) 



1.5 Beyond the Heston model 



17 



16.5 r 
16 

^ 15.5 

>. 

| 15 
o 

14.5 
14 
13.5 
13 



"D 
CD 

Q. 

E 




1 a 



10 25 



50 75 
Delta [%] 



90 



15.5 
15 

^ 14.5 



J? 
o 
> 

"a 
cd 

"a. 
E 



14 
13.5 

13 
12.5 

12 



-e — 1 M smile 
□ - Heston fit 



10 25 



50 75 
Delta [%] 



90 



17 r 

r-, 16 



= 15 

nj 

o 

I 14 

Q. 

E 

- 13 



12 L 



-e — 3M smile 
B - Heston fit 




10 25 



50 
Delta [%] 



75 90 



19 

18 

^ 17 
>. 

1 16 

o 
> 

"D 
CD 



I" 1 




10 25 



50 
Delta [%] 



75 90 



Figure 1.8: The EUR/USD market smile on July 22, 2010 and the fit obtained with the Heston 
model for r = 1 week (top left), 1 month (top right), 3 months (bottom left), and 6 
months (bottom right). 
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for some nonnegative deterministic parameter functions o~(t), K,(t), and 0(t). The formula for the 
mean turns out to be: 



E(«t) = g(t) = v e- K W + / n(s)6(s)e K ^- K M ds, 

Jo 

with K(t) = Jq k(s) ds. The result for the second moment is: 

E(v 2 ) = v 2 e~ 2K ^ + f {2K( S )0( S ) + a 2 (s)}g( S )e 2K ^ 2K ^ ds, 
Jo 

and hence for the variance (after some algebra): 

Var(« t ) = f a 2 (s)g(s)e 2K ^- 2K ^ ds. 
Jo 



(1.47) 



(1.48) 



(1.49) 
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Figure 1.9: The EUR/USD market smile on July 22, 2010 and the fit obtained with the Hcston 
model for r = 1 year (top left) and 2 years (top right). Again the term structure of 
the vol of vol and correlation visualizes the problem with fitting the smile in the short 
term (bottom panels). 
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The formula for the variance allows us to compute forward volatilities of variance explicitly. As- 
suming known values o^i and o~t 2 f° r times < T\ < T2, we want to determine the forward 
volatility of variance o-t x ,t 2 which matches the corresponding variances, i.e. 

4 2 ff(s)e 2K(s " T2) ds = (1.50) 

Tl 4 l5 (s) e 2K(s - T2) ds + r a 2 TuT2 g(s)e^- T ^ ds. 



The resulting forward volatility of variance is thus: 



? _ u\H(T % )-o\H(T x ) 
t t u t 2 H(T 2 )-H(Ti) ' 



(1.51) 
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where t 

H(t) = J g(s)e 2KS ds = |-e 2Kt + ±(v - &)e Kt + \{^~ «o) • (1-52) 

Assuming known values pt x and pT 2 for times < T\ < T2, we want to determine the forward 
correlation coefficient px t ,t 2 to be active between times T\ and T2 such that the covariance between 
the Brownian motions of the variance process and the exchange rate process agrees with the given 
values pt x and pt 2 - This problem has a simple answer, namely: 

Pt u t 2 = Pt 2 , T 1 <t<T 2 . (1.53) 
This can be seen by writing the Heston model in the form: 

dS t = St^dt + ^vidW^Y (1.54) 

dv t = k(6 - v t ) dt + pa^T t dW t {1) + y/l - p 2 a^dW t {2 \ (1.55) 

for a pair of independent Brownian motions and W^ 2 K Observe that choosing the forward 

correlation coefficient as stated does not conflict with the computed forward volatility. 

1.5.2 Jump-diffusion models 

While trying to calibrate short term smiles, the volatility of volatility often seems to explode along 
with the speed of mean reversion. This is a strong indication that the process 'wants' to jump, 
which of course it is not allowed to do. This observation, together with market crashes, has lead 
researchers to consider models with jumps (Gatheral, 2006; Martinez and Senge, 2002). Such 
models have been investigated already in the mid-seventies (Merton, 1976), long before the advent 
of SV. Jump-diffusion (JD) models are, in general, more challenging to handle numerically than 
SV models. Like the latter, they result in an incomplete market. But, whereas SV models can be 
made complete by the introduction of one (or a few) traded options, a JD model typically requires 
the existence of a continuum of options for the market to be complete. 

Bates (1996) and Bakshi, Cao, and Chen (1997) suggested using a combination of jumps and 
stochastic volatility. This approach allows for even a better fit to market data, but at the cost 
of a larger number of parameters to calibrate from the same market volatility smile. Andersen 
and Andreasen (2000) let the stock dynamics be described by a JD process with local volatility. 
This method combines ease of modeling steep, short-term volatility skews (jumps) and accurate 
fitting to quoted option prices (deterministic volatility function). Other alternative approaches 
utilize Levy processes (Cont and Tankov, 2003; Ebcrlcin, Kallsen, and Kristen, 2003) or mixing 
unconditional disturbances (Tompkins and D'Ecclesia, 2006). 
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